arXiv:1506.07062vl [cs.CV] 23 Jun 2015 


Improving Fiber Alignment in HARDI by 
Combining Contextual PDF Flow with 
Constrained Spherical Deconvolution 

J.M. Portegiesg R.H.J. Fick^, G.R. Sanguinetti^, 
S.P.L. Meesters^’^, G. Girard^’^, R. Duits^’^ 


1 Dept, of Math, and Computer Science, TU/e, The Netherlands 

2 Athena, INRIA Sophia Antipolis, Prance 

3 Academic Center for Epileptology Kempenhaeghe & Maastricht 
UMC+, The Netherlands 

4 SCIL, Universite de Sherbrooke, Canada 

5 Dept, of Biomedical Eng., TU/e, The Netherlands 
* E-mail: j.m.portegies@tue.nl 


Abstract 

We propose two strategies to improve the quality of tractography results com¬ 
puted from diffusion weighted magnetic resonance imaging (DW-MRI) data. 
Both methods are based on the same PDE framework, defined in the cou¬ 
pled space of positions and orientations, associated with a stochastic process 
describing the enhancement of elongated structures while preserving crossing 
structures. In the first method we use the enhancement PDE for contextual 
regularization of a fiber orientation distribution (FOD) that is obtained on indi¬ 
vidual voxels from high angular resolution diffusion imaging (HARDI) data via 
constrained spherical deconvolution (CSD). Thereby we improve the FOD as 
input for subsequent tractography. Secondly, we introduce the fiber to bundle 
coherence (FBC), a measure for quantihcation of fiber alignment. The FBC 
is computed from a tractography result using the same PDE framework and 
provides a criterion for removing the spurious fibers. We validate the proposed 
combination of CSD and enhancement on phantom data and on human data, 
acquired with different scanning protocols. On the phantom data we hnd that 
PDE enhancements improve both local metrics and global metrics of tracto¬ 
graphy results, compared to CSD without enhancements. On the human data 
we show that the enhancements allow for a better reconstruction of crossing hber 
bundles and they reduce the variability of the tractography output with respect 
to the acquisition parameters. Finally, we show that both the enhancement of 
the FODs and the use of the FBC measure on the tractography improve the 
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stability with respect to different stochastic realizations of probabilistic tracto- 
graphy. This is shown in a clinical application: the reconstruction of the optic 
radiation for epilepsy surgery planning. 


1 Introduction 


Diffusion weighted magnetic resonance imaging (DW-MRI) is a non-invasive 
technique for the characterization of biological tissue microstructure [^. In 
brain white matter, water molecules diffuse predominantly along axonal fibers. 
This results in an observable macroscopic orientation dependence in the DW 
signal, that is measured by scanning the tissue in multiple orientations and 
gradient strengths. To model the angular anistropy of the diffusion profile, 
diffusion tensor imaging (DTI) is widely used, but this has the limitation 
that only a single fiber direction can be estimated per voxel [^ . It is estimated 
in that more complex fiber configurations occur in approximately 90% of the 
white matter voxels. To overcome this, high angular resolution diffusion imaging 
(HARDI) techniques are used, that can describe more complex (crossing) fiber 
configurations. An overview of HARDI techniques can be found in [^. Here we 
use the method of constrained spherical deconvolution (CSD) [^, that from the 
initial diffusion data constructs a fiber orientation distribution (FOD), which 
models the distribution of fibers along different directions. 

Tractography methods are often used in the DW-MRI pipeline to provide in¬ 
sight in the structural connectivity of the white matter bundles. Independently 
of the model used for interpreting the DW-MRI data, noise originating from the 
scanner, acquisition artifacts and partial volume effects are likely to result 
in spurious (aberrant) fibers in the tractography output. To improve the data 
on which the tractography is performed, different regularization methods can 
be used. Methods exist that apply filtering for the reduction of noise directly 
on the DW-MRI data [8 - 10 , other methods aim to regularize the DTI tensor 


fields 11 ■ 15 


ual voxels 16 -18 


On HARDI data the regularization can be performed on individ- 
or in combination with the local spatial information 
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We introduce two new strategies based on the same underlying principle to 
improve fiber alignment in tractography results, in order to have more reliable 
information on the structural connectivity of brain. First we perform contex¬ 
tual regularization to the FOD obtained with CSD, see Fig. and secondly 
we introduce a fiber to bundle coherence (FBC) measure that can be applied 
to any fiber bundle to classify and remove spurious fibers, see Fig. 03- Both 
approaches are based on a partial differential equation (PDF) framework intro¬ 
duced in 24 ■ 27 , where the Fokker-Planck equation of a stochastic process for 
enhancement of elongated structures is considered. These type of PDF-based 
enhancement methods have been widely used for the processing of 2D-images. 
In this framework, images are represented in the extended space of positions 
and orientations via a stable invertible orientation score 28 , that associates 


to every location an orientation distribution of the local image features (lines 
and contours). Then, the stochastic processes for contour completion [^ 34 
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Figure 1: The proposed pipeline of the paper. CSD is used to estimate 
an FOD from DW-MRI data. The FOD is enhanced (A) with PDF techniques. 
Then a deterministic or probabilistic tractography is applied to the (enhanced) 
FOD (probabilistic shown here, with coloring indicating the fiber direction). 
In the lower right figure, we applied our coherence quantification method (B), 
based on the same PDF framework, which shows that blue fibers are well aligned 
(high Fiber to Bundle Coherence (FBC)) and yellow fibers are spurious (low 
FBC). The green arrows indicate the steps in which the contextual PDFs are 
used. 


and contour enhancement 35-37 (see Fig. [^) on this extended space xi S' 
induce crossing preserving completion and enhancement of lines 28 


The DW-MRI data that we use is naturally defined on the coupled space 
X S^ of 3D positions and orientations. As in the 2D case, crossing preserving 
enhancement of line structures is required, for which we use the 3D extension 
of the 2D stochastic process for contour enhancement, introduced in 25 . The 


linear PDF corresponding to this stochastic process can be solved by convolution 
of the initial condition with the kernel of the PDF. This kernel is also a function 
on the position-orientation space and can be seen as a transition distribution 
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Figure 2: Stochastic interpretation of the contour enhancement ker¬ 
nels. A. Accumulation of 300 sample paths drawn from the underlying stochas¬ 
tic process of the contour enhancement PDE in x 5^, projected on the xy- 
plane. B. The contour enhancement kernel arises from the accumulation of 
infinitely many sample paths. The gray-scale contours indicate the marginal 
of the kernel, obtained by integration over 5'^, the red glyphs are polar graphs 
representing the kernel at each grid point. C. The contour enhancement kernel 
oriented in the positive z-direction in x S'^ can be visualized on a grid with 
glyphs that in this case are spherical graphs. 


from the origin (in position and orientation) to neighboring elements. From 
the stochastic point of view, the kernels arise as limits of the accumulation of 
infinitely many sample paths drawn from the stochastic process, illustrated in 
Fig. [^. For mathematical details of the underlying stochastic processes of the 
PDFs, see §10.1]. The general idea needed for this article is sketched in Fig. 

In Figs. and we show the contour enhancement kernel using glyph 
visualization on a grid, each glyph being a polar (red, 2D) or spherical (blue, 
3D) graph plot where in every orientation the (spherical) radius is proportional 
to the value of the kernel. This type of visualization is used throughout the 
paper for functions defined on the space of positions and orientations. 

Recently, many authors [23p5p^34|38f|4^ demonstrated the advantages of 
contextual processing of DW-MRI data. The general rationale behind contex- 
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tual processing is to include alignment of local orientations and their surround¬ 
ings (i.e. the context) on the coupled space of positions and orientations. For 
this alignment of local orientations, roto-translations are needed, which imposes 
a non-Euclidean structure in the PDE-based processing as we explain in Section 
|2.2| More details on the embedding of x 3“^ in the roto-translation group 
S'if(3) can be found in 25 . This demonstrates how either the completion or 


enhancement PDEs can be used to extrapolate DTI information to increase the 
angular resolution and resolve some fiber crossings. This idea was shown to be 
promising in clinical experiments 38||^ , but in some cases extreme parameters 
had to be set to obtain clear maxima at crossings (where DTI data is inade¬ 
quate). Therefore in this paper we introduce and test the combination of CSD 
with contextual enhancements. The method proposed in [34| uses an advection- 
diffusion equation (that we called contour completion above) to improve HARDI 
data to obtain connectivity measures. In our work we rely on a purely diffusive 
process, contour enhancement, which in contrast to contour completion does not 
suffer from singularities 27 and is less sensitive to small perturbations of the 


initial conditions. This property makes the enhancement process more suited 
to be combined with the sharp angular distributions produced by CSD. As the 
methods mentioned above still result in broad angular distributions, they need 
to be combined with some sharpening method. To this end, a geometric mor¬ 
phological sharpening based on erosions was presented in 23 27 42 . Another 


related method presented in 40,41 is the so-called fiber continuity model in 


which purely spatial regularization is considered in combination with spherical 
deconvolution as alternative to the non-negativity constraint in the classical 
CSD 43 . In Section 2.2 we demonstrate the importance of including also an 


angular regularization term. 


1.1 Contributions 

The first contribution of this article is to study the combination of the widely 
used CSD method with a regularization induced by the enhancement PDE act¬ 
ing on the FOD. Since the FOD obtained with CSD consists of sharp angular 
profiles, it is well-suited as an initial condition for the enhancement PDE, that 
typically has a smoothing effect on the orientation distributions. The contex¬ 
tual regularization method reduces non-aligned crossings in the FOD, allowing 
for a better alignment of fibers when tracking is applied on the enhanced FOD. 
We show that this method is therefore useful to reduce the number of false 
positive fibers, but mainly to find more true positives in the tractography out¬ 
put. Although in this paper we compare to the classical CSD method, the PDE 


enhancements can also be applied to extensions of this method 44 - 48 


The second contribution of this article is to introduce the fiber to bundle 
coherence (FBC) measure. The motivation for this measure is that, especially 
probabilistic, tracking methods typically produce spurious fibers that should be 
removed from the tractography. In contrast to the first approach, this method 
serves as a post-processing tool. For the computation of the FBC we regard the 
fiber bundle as a set of oriented points, by considering for every fiber point also 
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the local tangent to the fiber. We construct a density using the enhancement 
PDE with an initial condition that is a sum of superposed (^-distributions at 
every oriented point in the bundle. The construction of such a density from 


tracks relates to track density imaging 49 and track orientation density imaging 


50 , though here the use of the contour enhancement kernels, Fig. allows to 
use a sparse set of fiber tracks. The FBC, a measure for spuriousness of fibers, 
is computed by efficient integration of this fiber-based density. Fibers that are 
most spurious according to the FBC can be removed from the tractography, 
resulting in a better aligned fiber bundle. Complementary to the first method, 
this FBC measure has the purpose to remove false positives in a tractography. 


1.2 Structure of the Article 

Section covers theory of the individual parts of the pipeline as outlined in 
Fig. [TJ consisting of CSD, PDE enhancements, tractography and coherence 
quantification in Sections 2.1||2.4[ respectively. In Sectionj^we provide extensive 
validation of the combination of CSD and PDE enhancements and the FBC, 
using three experiments: 


1. First we use the Tractometer evaluation system 51 52 on the ISBI 2013 
HARDI reconstruction challenge dataset 53 , a digital phantom with 


known ground truth, to demonstrate how contour enhancement improves 
both the local FOD reconstruction and the global connectivity of fiber 
bundles compared to CSD, see Section [3T| 


2. In Section 3.2 we show on a human DW-MRI dataset, containing different 
crossing bundles, that CSD combined with enhancements yields an FOD 
that is more robust with respect to the 6-value and the number of gradient 
directions used in the acquisition. Furthermore, we make a comparison 
with earlier work involving erosions and nonlinear diffusion of FODs di¬ 
rectly applied to a DTI-model 23 27 , that was based on the same data. 


We show that with our method the glyphs are sharper at the locations 
where bundles cross. 


3. Finally in Section [3.3[ we show an experiment with clinical data in which 
we reconstruct the optic radiation (OR) to determine the position of the 


tip of the Meyer’s loop, that is of interest in epilepsy surgery planning 23 


54 -57 . Accurate estimation of this position is difficult due to the presence 


of spurious fibers in the reconstruction of the OR. We show that both the 
FOD enhancement and the FBC measure, see Fig. and in particular the 
combination of the two allow for a more stable determination of the tip of 
the Meyer’s loop. Here ‘more stable’ means less variation with respect to 
stochastic realizations in the probabilistic tractography results. 


Conclusions and a discussion can be found in Section |4| 
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2 Methods 


In this paper it is assumed that we have HARDI data as input, from which 
we derive an FOD U that models the orientation of fibers in each voxel, i.e. 
U X R+. For this we use CSD concisely described in Section 

|2.1[ as it gives sharp angular profiles and is able to distinguish multiple fiber 
directions within a voxel. 

Then we use the enhancement PDF for diffusion of the FOD U, coupling spa¬ 
tial and angular information. The combination of CSD and such enhancement is 
a powerful method to obtain an enhanced FOD in which the coherence inherent 
in the data is included, providing a more coherent input for the tractography. 
The enhancement technique is explained in Section |2.2[ 

We use the MRtrix algorithm for both deterministic and probabilistic 
tractography to estimate the structural connectivity in the brain. In the de¬ 
terministic tractography, fiber tracks are obtained by integrating a directional 
field, given an initial position and direction. The directional field is given by 
the locally maximal orientations in the glyphs. In contrast to deterministic 
tractography, the probabilistic tractography method of MRtrix samples the ori¬ 
entations from the entire FOD and does not use just the maxima. More difficult 
paths can be reconstructed than with deterministic tracking, but typically also 
many spurious fibers are produced due to the probabilistic sampling. Both 
the deterministic and the probabilistic method are explained in more detail in 
Section 12.31 


In Section 2.4 we introduce our new technique to quantify the coherence of 
fibers with respect to all the fibers in a bundle, based on the same PDF theory 
as employed for the contextual enhancement in Section |2.2[ We explain how 
the kernel of the enhancement PDF is used to construct a tractography-based 
density, how the FBC is computed and how this measure is able to classify 
spurious fibers in a tractography. 


2.1 A Brief Review of CSD 


In CSD it is assumed that at each voxel position y the measured signal Sy : 

—)■ K can be represented by a spherical convolution of the FOD /y : —>■ K 

with a response function AT : —)■ M, that is estimated from the data 58 . Since 


the spherical deconvolution to determine the FOD is ill-posed, a non-negativity 
constraint is included as in [4^[^. Then, given the signal S'y(n) for a sample 


of orientations n G the solution of CSD is found by iteratively solving the 
minimization problem: 


/^+^(n) = argmin [ \{K g){n) - Sy{n)\^da{n) +\^ [ |(£j. (g))(n)pdcr(n), 

ogL2(S2)./s2 Js^ 


■v* 

Data Driven 


Regularization 


( 1 ) 
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for i = 1,... ,imax) with imax the maximum number of iterations. Here K £ 
L 2 ( 5 '^) is aligned with and symmetric around the z-axis, the convolution *52 
is the usual spherical convolution 
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da{n) is the Jacobian of the sur¬ 
face measure in orientation n and A is a parameter to influence the trade-off 
between the data driven term and regularization term. The linear operator 
Ch ■ L 2 (S'^) —>■ L 2 ( 5 '^) in the regularization term gives the non-negativity con¬ 
straint and is defined by: 


{Chf){n) = f{n)H{Th - h(n)), for given h G L 2 (S'^), 


( 2 ) 


where H is the Heaviside function and Th is a threshold equal to a fixed factor 
r times the mean of h. The initial function /y for the iteration is computed 
by taking only the data driven term of Eq. Q. The iteration stops when 
successive iterations yield the same result, typically after 5 to 10 iterations 
Throughout the paper, we call U the FOD obtained by 
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C^(y,n) = /y““(n). 


(3) 


In practice CSD is performed using spherical harmonics with a maximal 
spherical harmonic order of 8 (/max = 8 ) as discussed in 61 . 


Improvements to the original CSD exist to modify and improve the response 
function, either by recursive calibration or auto-calibration 
multiple acquisition shells 


44 45 


46 or by including anatomical data 47 48 


by using 
The 


latter two methods aim to reduce the partial volume effects, where CSD is likely 
to produce spurious fiber orientations. These partial volume effects can occur 
when in a voxel multiple tissues or multiple bundles with different orientation 
are present. Here we use the classical CSD as it is the basic technique available 
in several neuroimaging packages. However, we stress that our method is not 
restricted to this type of CSD. In any case, our method aims to reduce non- 
aligned crossings in the FOD, also the ones induced by partial volume effects, 
as we will show in several experiments in this paper. Further improvement of 
the methodology can be expected when including recently extended and more 
elaborate CSD techniques 45 46 48 , but this is left for future work. 


2.2 Contour Enhancement (Step A) 

To improve alignment of neighboring glyphs of the FOD U, recall the glyph field 
visualization in Figs. 0 and[^, we apply contextual enhancements. Before we 
specify the PDF we consider for this enhancement, we first need to express the 
notion of alignment in mathematical terms. To this end, let us consider Fig. 
where it is shown that the notion of alignment cannot be supported by a decou¬ 
pled, flat Cartesian product x with the combined Euclidean distance. It is 
clear that the green bar at (yi, rii) is better aligned with the gray bar at (yo, no) 
than the orange bar at (y 2 , n 2 ), even though the distances in the space x 

are equal, i.e. yc/R 3 (yo,yi) +c/| 2 (no,ni) = (yo> yz) + (>^o, n 2 ). This 

means that in order to appropriately describe the concept of alignment, we must 




























consider more than just the amount of spatial displacement and the amount of 
change in orientation. Coupling these two types of motion (via rigid body mo¬ 
tions) is a solution to this problem 25 . The coupling follows very naturally by 


expressing the motion of an oriented particle (y, n) in terms of a moving frame 
of reference determined by its orientation. That is, spatial movement along 
the orientation n should be much cheaper than spatial movement in the plane 
orthogonal to n. This creates a natural anisotropy for spatial movement. For 
angular motion we need isotropy. This extra structure can be obtained by em¬ 
bedding the space of positions and orientations in the rigid body motion group. 
This means that an element (y, n) € x is identified with the rigid body 
motion (y,Rn), where Rn is any rotation matrix such that Rne^ = n, with 
Bz G S'^ pointing to the north pole. We denote this space of coupled positions 
and orientations by x so we have 


x52 9(y,n)o(y,R„)GK3^52 


( 4 ) 


The group x is equipped with the following (non-commutative) group 
product: 


(y,R)(y',R') = (y + Ry',RR')- 


( 5 ) 


This product moves oriented elements in a shift-twist fashion, rather than by a 
rotation followed by an independent translation. Due to this shift-twist group 
product in Eq. ([^, we automatically express motion of oriented particles in 
terms of a moving frame in x which makes this space well-suited for the 
application of our contextual enhancements. Nevertheless, in the remainder of 
this article this space can still be regarded as the Cartesian 5D space x 
where we secured the coupling of positions and orientations via our specific 
choice of differential operators and diffusions that are applied. 

To improve alignment of FOD glyphs, we use a particular diffusion process 
called contour enhancement that uses both spatial and angular diffusion in the 
extended space of positions and orientations 25 . Given a structure (think of 


a fiber bundle) in this space, see Fig. we apply spatial diffusion only in the 
direction of the structure, not in the spatial plane perpendicular to it. Angular 
diffusion is applied in the plane tangent to 5^ at the point n. This diffusion 
process enhances elongated structures, while preserving crossing structures, and 
is given by a Fokker-Planck type of system, a linear diffusion equation on x 5*^. 
For t > 0, y G n G this system can be expressed as: 


f dtW{y,n.,t) = {Dss^n ■ 
\ 4F(y,n,0) = U{y,n). 


+ D 44 AS 2 ) W{y,n,t), 


( 6 ) 


Here W (y, n, t) is a scale space representation in (K^ x S'^) xK+ 62 . The symbol 


Vy denotes the gradient with respect to the spatial variables and A 52 is the 
Laplace-Beltrami operator on the sphere. Parameters D 33 > 0 and D 44 > 0 are 
related to the amount of spatial and angular diffusion, respectivly. Parameter 
t > 0 is the diffusion time of the contour enhancement process. It can be seen 
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(yi.ni) 



(y2,n2) 


Figure 3: The concept of alignment requires a coupling of positions 
and orientations. The pair of position and orientation (yo,no) is better 
aligned with (yi,ni) than with (y 2 ,n 2 ), even though spatial and angular dis¬ 
tances are equal. Formally we can say that the sub-Riemannian distance on 
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is smaller between (yo,no) and (yi,ni). 




Figure 4: Diffusion in the space of positions and orientations. Spatial 
diffusion is applied in the direction of the fiber (left), angular diffusion is applied 
in the tangent plane perpendicular to the fiber direction (right). 


as a Brownian motion process, recall Fig. where particles are allowed to 
spatially moye back and forth in the direction they are heading, or change their 
direction, but are not allowed to step aside (comparable to the moyement of a 
car). 

We refer to the solution of Eq. (|^ as the enhanced FOD. It can be obtained 
via a finite difference scheme [^, or via a convolution with a kernel pt : x 

5*2 ^ M+: 
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(7) 


W{y,n,t) = {pt *h3x,s2 U) (y,n) 

= f f Pt((Rn')^(y - y'))R-n'n) • t/(y',n') dy'dcr(n'). 

Js^ JR3 

A basic approximation to the exact Green’s function of the contour enhancement 


PDE is known 


25 


and can be written as the product of Green’s functions 




in the following way: 


Pt(y,n) = ■^D 33 t^/J^il^■pf'^^^\z/ 2 ,x,| 3 )■pf^^\z| 2 ,-y,y), (8) 


with n = n(/3,7) = Re^, 7 Rej,,^e^ = (sin /3, — cos /3 sin 7 , cos /3 cos 7 )^, /3 G 
[— 7 r, 7 r), 7 G (—f, §)• The xi kernels are given by 


with 


pf ^^\x,y,9) = 


1 


327 rt ^ 044^1)33 




EN{x,y,e) 


(9) 


EN{x,y,e) = 



0/2 


To avoid numerical errors, we use the estimate 


e /2 


tan( 0 / 2 ) 
( 10 ) 

tan(e/2) ~ l-(i2/24) 1^1 < fs' 

This approximation is easy to use and allows for efficient implementation [64| . 

From the approximation kernel in Eq. Q it can be seen that problems could 
occur when either D 33 = 0 or D 44 = 0. To this end, a necessary and sufficient 
condition for the existence of a smooth solution kernel for the evolution process 
in Eq. § is given by the Hdrmander requirement 
more general situations than the one here, see e.g. 


y 


65 
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This condition applies to 
but for the specific case 
of contour enhancement the requirement is satisfied iff 033,044 > 0. Setting 
I ?44 = 0 would result in a singular non-smooth kernel, which has numerical 
disadvantages. More importantly, apart from this theoretical issue the need for 
both spatial and angular diffusion can also be argued from a practical point of 
view, as is illustrated in Fig. We use an artificial example in which a curved 
fiber bundle is present, shown in the left figure. When the input is diffused with 
O 44 = 0 as in the middle of Fig. the peaks stay distinct and point in the 
wrong direction. On the other hand, when O 44 > 0 as in the right figure, due 
to the angular diffusion the peak is redirected and the glyphs lie better aligned 
with the fiber bundle. Hence O 44 > 0 is needed to ensure the crucial interaction 
between different orientations. Finally we recall the relation between Tikhonov 
regularization and diffusion, see e.g. 25 Thm 2], which allows us to connect 
diffusion with £>44 = 0 with the fiber continuity model in 40 41 . This model 


does not suffer from the inconvenience of considering only spatial regularization. 
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Figure 5: The importance of including angular diffusion. (Left) artificial 
input data of a curved bundle. (Middle) diffusion with II 33 > 0, ZI 44 = 0. 
(Right) contour enhancement with 1133,1144 > 0. Fiber propagation with ZI 44 = 
0 leads to crossing artefacts rather than smooth fiber enhancement. 


as they represent the FOD in a truncated spherical harmonic basis. When the 
enhancements are used in combination with probabilistic tractography, we first 
apply a standard sharpening deconvolution transform to the FOD as described 
to maintain the sharpness of the FOD. 
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2.3 Tractography 

As the next step in the pipeline we use the MRtrix tractography algorithm , as 
implemented in http://www.brain.org.au/software/index.html^mrtrix, version 
0.2.12. It allows us to perform deterministic and probabilistic fiber tracking 
on spherical harmonic representations of the (enhanced) FOD. To have a fair 
comparison between trackings on the FOD and the enhanced FOD, we use the 
parameter settings as explained next. 

• In the deterministic tracking of MRtrix, seed points are randomly selected 
from a seed region. The initial direction is sampled randomly and every 
next step follows the direction of the most aligned FOD maximum. If this 
maximum is below a threshold value, the fiber terminates. This threshold 
(cutoff) is set to 10% of the maximal angular response of the FOD. There 
is no constraint on the maximal curvature of the fibers. To prevent that 
fibers have an initial direction that is not aligned with the fiber bundle, 
we force the initial direction to be approximately in the direction of the 
maximal FOD peak, by setting the initial cutoff to 0.9. The step size is 
set to 1/lOth of the voxel size as is suggested in [^. Tracks proceed in 
both directions from the seed point and terminate either when they hit the 
boundary of the volume or mask (if applicable), or due to the threshold 
stopping criterion. 

• In the probabilistic case, starting from the seed region, every next step 
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follows a direction randomly sampled from the FOD. Here we set the 
minimal radius of curvature to 1 mm, the default value in the MRtrix 
algorithm. Optionally, a target region of interest is used to select only 
those fibers that cross this region. 

We base our choice for deterministic tractography or probabilistic tracto- 
graphy on the application. If only a seed region is specified, as in Sections |3.1| 
and |3.2[ we use deterministic tractography. In this case there is too much free¬ 
dom in the probabilistic algorithm and the streamlines show a lot of spurious 
behavior. Here, a probabilistic approach could make sense if extreme amounts 
of tracks are used for track density methods. As we do not pursue these meth¬ 
ods here, we prefer to use deterministic tractography. If both a seed region 
and an end region are specified, as in the optic radiation application in Section 
|3.3[ we prefer to use probabilistic tractography. It is known that deterministic 
tractography in this case provides only a few of the possible pathways from the 
seed to the end region, whereas reconstructions with probabilistic tractography 
are much fuller. 

Probabilistic tractography results typically contain many false positive fibers. 
Streamlines that are anatomically implausible can be removed with scoring 
methods [23[|56| or by imposing anatomical constraints. Even when using these 
methods, the filtered tractography output can still contain fibers that deviate 
from the fiber bundle and are likely to be spurious. In the next section, we pro¬ 
pose a coherence measure for fibers in a fiber bundle in order to classify these 
spurious fibers. 

2.4 Coherence Quantification of Fiber Bundles (Step B) 

In this section we introduce our second contribution of the paper, a fiber to 
bundle coherence (FBC) measure to quantify the coherence of each fiber with 
respect to all other fibers in the bundle, recall Fig. [^. A spurious fiber, as 
schematically shown in Fig. is isolated from or poorly aligned with the bulk 
of the tracks and is therefore unlikely to represent the underlying brain structure. 
Fibers with low coherence, i.e. a low FBC, can then be classified as spurious. 

To classify a fiber as spurious, we first construct a density by regarding 
each fiber as a superposition of ^-distributions in xi S'^ and convolving this 
distribution with the kernel in Eq. Q. This density is independent of the 
underlying data and is based purely on the collection of fibers F. Integration 
of this density along a part of length a of a fiber gives a local measure for the 
coherence of that part. 

Next we explain the mathematical techniques that support the idea in Fig. 

We denote the fibers from a tractography output by yi(s) G K^, 1 < j < N, 
0 < s < li, with s the arc length parameter, k the total length of fiber i and N 
the number of fibers. Now let ni(s) = yi(s) be the tangent of the fiber, so that 
7i(s) = (yi(s), ni(s)) forms a curve (fiber) in x S'^. By construction, ni(s) 
points in the forward direction of the fiber. Since in DW-MRI data antipodal 
orientations are identified, we also consider 7i(s) = (yi(s), —ni(s)). The com- 
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Figure 6 : Construction of the LFBC. The local fiber to bundle coherence 
(LFBC) is constructed for a set of fibers (gray lines), illustrated in 2D for 
simplicity, as follows. Every local tangent in the tractography contributes to 
the density, by considering it as a (^-distribution. We convolve this with the 
contour enhancement kernel, shown on the left for two points, visualized as 
in Fig. and with the coloring indicating the contribution of the kernel to 
the LFBC. Doing so for all points, fiber points that are isolated from or badly 
aligned with other fibers receive low contributions, such as the outlying fiber. 
The LFBC along the fibers is displayed on the right. 


plete fiber bundle is defined as T := {yi | i = 1,..., N} U { 7 ^ | i = 1,..., N}. A 
discrete formulation of a fiber i with Ni points is given by: 


if :=7i(s^) = (y*(s^),nz(s^)) =: (y-,n^), 


j - 1 
N,-l 




( 11 ) 

This way there are Ntot = 2 X)i=i elements in T. Now we regard every 
point as a (5-distribution in x centered around A density for 

the entire bundle is then constructed as follows: 


-Fr(y,n) 


1 

Aftot 


2 N Ni 


a—li—1j—1 


( 12 ) 


with index j running over points within a fiber, i running over all fibers and a 
taking care of including forward and backward orientations. We use the same 
evolution process as in Eq. (§ in which F = now serves as initial condition, 
to create a diffused density (y, n) (->• Wpiy, n, t): 


f dtWpiy, n, t) = (D 33 (n • Vy)^ -k D 44 AS 2 ) Wpiy, n, t), 

I WF(y,n, 0 ) = F(y,n). 

We solve the system in (13) by convolution with the corresponding kernel, recall 
Fig. and call this the local FBC (LFBC): 


LFBC(-,r) = Wpi-it) = {pt F)(.), (14) 
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with the shift-twist convolution as given in Eq. Q. This is illustrated in Fig. 
[^in the 2D case. We can now define the FBC for fiber 7 ^ with respect to the 
bundle T as the integral of this density along the fiber: 

FBC( 7 „r) = -/ LFBC( 7 ,(s),r)ds. (15) 

k Jo 

This results in a global property of the fiber, but spurious fibers often only 
locally deviate from the bundle as in Fig. To this end, we compute for each 
fiber the minimum of such integrals along the fiber over intervals of length a: 

2 pa-\-OL 

FBC“( 7 i,r) = min — / LFBC( 7 i(s),T) ds. (16) 

£i6[0,7 —q] O . 

The parameter a defines the scale over which spuriousness of fibers can be 
detected and is much smaller than the average fiber length. Our primary interest 
is not the FBC“ value itself, but rather how it compares to the average coherence 
of fibers in the bundle, so finally we define the relative fiber to bundle coherence 
(RFBC) as: 


RFBC( 7 „r) 


FBC“( 7 „r) 
AFBC(r) ■ 


(17) 


Here AFBC(r) is the average fiber to bundle coherence indicating the overall 
coherence of the N fibers in the bundle F, defined as 


N 


AFBC(F) = -5]FBC(7.,r). 


(18) 


i=l 


To summarize, the RFBC( 7 i, F) of a fiber 7 ^ in a bundle F is a measure for how 
well aligned the least aligned part of % is, compared to the average coherence 
of the total bundle. 


In practice, we evaluate the convolution in Eq. (14) only in the fiber points. 
We compute the LFBC( 7 (‘,F), the diffused density in the oriented point = 


(y"-',!!*), recall the notation in ( 11 ), as follows: 


LFBC( 7 f,F) 


N Ni 








cr=l j = l q=l 


(19) 


where R„i is any rotation matrix such that R ; = n(, index q sums the 

contributions along a fiber, index j runs over all the fibers and a as before. The 
FBC“ can then be computed as follows: 


, a+a 

FBC“( 7 .,r)= min - ^ LFBC( 7 f,F), 

aG\0.Ni—Oi] O. 


( 20 ) 


15 






where a, a G N in this discrete case, so the LFBC is summed along short intervals 
of the fiber. Likewise, the AFBC can be computed as 

W TV, 

AFBC(r) = -^ —^LFBC( 7 f,r). (21) 

i=l ^ * fc=l 

We apply this method in Section |3.3| for quantifying the coherence of tracto- 
graphy results of the optic radiation and classifying the spurious fibers. 


3 Experiments and Results 


In this section we extensively test the performance of our CSD enhancement 
method (A) and the FBC method (B), recall Fig. [^and Sections 2.2 and 2.4 
in three different experiments: 


We use the HARDI Reconstruction Challenge dataset 66 , which is arti¬ 


ficial data with known ground truth, to quantitatively evaluate the CSD 


enhancement method (A) on deterministic tractography in Section 3.1 


In Section 13.21 we show on DW-MRI human brain data that the enhance¬ 
ment (A) have a positive effect on deterministic tractography, for different 
acquisition protocols of the data. Furthermore, on this DW-MRI dataset 
and on the phantom dataset we compare our method to previous work 27 


where a DTI-based FOD is used in combination with nonlinear PDF flow. 

• In the third and last experiment, we reconstruct the optic radiation in 

We include an extensive evaluation 


human clinical data, see Section 3.3 


of our methods, the enhancement of the FOD (A) and the use of the FBC 
to classify and remove spurious fibers (B), and the combination of both 
methods. We show that the reproducibility of the probabilistic tracto¬ 
graphy has increased, resulting in a more stable localization of the tip of 
the Meyer’s loop. 


For all datasets Mathematica 67 was used to perform the contour enhance¬ 


ment algorithm and the CSD, which in practice produces the same results as the 
MRtrix CSD implementation when the same deconvolution kernel is used. MR- 
trix software was used to perform fiber tractography. The coherence quantifi¬ 
cation was implemented in C-I--I-. In Section [3.1| we make use of the Tractome- 
(http://www.tractometer.org/) to evaluate tractography results. 


51 52 


ter 

Visualization was done in either the FiberNavigator (https://github.com/ 
scilus/f ibernavigator, 68 ), Mathematica, or the open source vIST/e tool 


(Eindhoven University of Technology, Imaging Science & Technology Group, 
http://bmia.bmt.tue.nl/software/viste/). 
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3.1 HARDI Reconstruction Challenge 


The following experiment is performed on a digital phantom dataset that was 
designed for the ISBI 2013 Reconstruction Challenge 


combination with the Tractometer 


66 69 . It is used in 


as a benchmark to compare different 


reconstruction and tracking methods. The phantom is inspired by the Numerical 
Fiber Generator 70 and the code to reproduce it is freely available as part of the 
Python package Phantomas (http://www.emmcLnuelcaruyer.com/phcLntomas. 
php). This synthetic dataset is of size 50 x 50 x 50 voxels with a resolution 
of 1 X 1 X 1 mm^. It consists of 27 simulated white matter bundles, designed 
to resemble challenging branching, kissing and crossing structures at angles 
between 30 and 90 degrees, with various curvature and bundle diameters ranging 
from 2 mm to 6 mm. An image indicating the ground truth fiber configuration 
is shown in the centre of Fig. 

The idea behind the signal simulation is that every voxel is subdivided into 
multiple sub-voxels, each one with its own attenuation profile. The final signal 
arrives from integrating the contribution of all the sub-voxels. Then, it is pos¬ 
sible to combine multiple compartment types in every voxel with added Rician 
noise. This allows for modelling complex configurations as well as taking into 
account partial volume effects. While the Numerical Fiber Generator uses a 
tensor-like model to simulate the signal in the sub-voxels, Phantomas uses a 
GHARMED-based mode l [tI] . The GHARMED model based on the Soderman- 


Jonsson cylinder model 72 captures well the non-Gaussian behaviour of the 


diffusion signal for large b-values. The main reason why we selected the ISBI 
phantom is that it is linked with the Tractometer that allows for performing 
quantitative evaluations of the tractography results, using global metrics as 
demonstrated in the subsequent experiments. 

For the experiments presented in this section we used 64 uniformly dis¬ 
tributed gradient directions using a 6-value of 3000 s/mm^ with different signal 
to noise ratios (SNRs). We use spherical harmonics in GSD with maximal order 
8 , resulting in 45 estimated coefficients on each position. We then enhance the 
resulting FOD functions using our contour enhancement algorithm with varying 
parameters. From the evolutions described in Eq. (|^ we see by a basic rescaling 
argument that it is sufficient to vary t and the ratio 1 / 33 / 1144 . The larger this 
ratio, the more preference the spatial diffusion gets over the angular diffusion, re¬ 
sulting in elongated kernels (visualized by thin glyphs). A smaller ratio I/33/7/44 
is better suited in regions where the curvature of bundles is higher (visualized by 
thicker glyphs). The higher the diffusion time t, the more context is taken into 
account. When t is too large, fiber bundles with high curvature can be damaged 
or false positives could be created. Taking this into consideration, we choose our 
parameters as follows: we fix spatial diffusivity parameter I/33 = 1.0, we take 
the angular diffusivity parameter Z/44 e {0.005,0.01,0.02,0.04} and diffusion 
times t G [0, 5]. 

Tractography results for the entire dataset are shown in Fig. We can 
recognize the positive effect of the enhancements on deterministic tractography: 
we see less dropouts, better aligned fibers and better continuation of fibers 
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Enhanced CSD 



Figure 7: Tractrography results on the ISBI Challenge dataset. De¬ 
terministic tractography results of CSD (left) and enhanced CSD (right) with 
SNR = 2 (top) and SNR = 4. The colors correspond to the direction of the 
fibers. The dataset consists of crossing, branching and kissing fiber bundles. 
The tractography on enhanced CSD results in better aligned hbers and a fuller 
reconstruction of the bundles. The ground truth configuration of the bundles is 
depicted in the center. 


at crossings. An extensive quantification of the performance of our method 
is done at the voxel level using the FODs and at the macroscopic level using 
tractography in Sections [3.1.1| and |3.1.2[ respectively. Both sections support the 
results summarized in Fig. 

3.1.1 Local Metrics 

We compare reconstructed FODs locally with the ground truth using only the 
orientation of the peaks. Let M be the set of voxels v in the white matter mask, 
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then we denote the ground truth number of peaks in a voxel v by N'" and the 
orientations corresponding to the peaks by i = 1 ,..., N^. 

Maxima of the constructed FOD are found by evaluating the FODs on a 60*^ 
order icosahedron tessellation with 18606 antipodally symmetric points, giving 
an angular resolution of less than 1 degree. Maxima are taken into account 
only if it exceeds a threshold of 0 . 1 , the same value we use as threshold in the 
tractography. Let be the set of peak orientations in voxel v estimated from 
the FOD. The average angular error in degrees can then be computed by: 


d = 


AT” 

E E min 

veM i=l 


180 


arccos n; 


E 

vGM 



( 22 ) 


In the top row of Fig. [^we show the effects contour enhancement for different 
ratios of ZI 33 and ZI 44 upon variation of the diffusion time. The results are 
given for substantially low SNR levels 10, 6 and 4 and 2. These SNRs are 
computed w.r.t. the non-DW image. Specifically, if the b=0 intensity is 1 then 
the standard deviation of the Rician noise distribution is 1/SNR. In all cases a 
clear improvement is found compared to CSD without enhancements and the 
more noise, the more the angular error is decreased. Higher diffusion times give 
better results and around t = 5 the angular error is almost stable. It can also 
be seen that the combination of CSD with enhancements at lower SNRs gives 
lower angular errors than just CSD for the higher SNRs. 

There is no significant difference in the FODs between the different D 44 val¬ 
ues. Even though it is visible that more angular diffusion leads to fatter glyphs, 
for the orientation of the peaks the precise value of D 44 is not of great impor¬ 
tance: the angular errors for D 44 = 0.005 are slightly smaller, but there is not 
much difference with the higher values of D 44 . 


3.1.2 Global Metrics 


At the macroscopic level we are interested in the impact of the enhanced lo¬ 
cal reconstruction on the quality of the global connectivity. The deterministic 
MRtrix tractography is used as described in Section |2.3[ with seeds randomly 
selected in the white matter mask. The tracks have a minimum length of 10 
mm and new seed points are chosen until 10000 streamlines are selected. For 
every FOD, the tractography is repeated five times with the same settings, to 
average out the variability in the tracking algorithm output. We then use the 


Tractometer 52 to perform a fiber tracking analysis based on the ground truth 


and the five results are averaged. The Tractometer outputs values for various 
metrics, from which we use the Valid Connections (VC), Invalid Connections 
(IC) and No Connections (NC). They indicate the percentage of tracks that 
correctly connect, incorrectly connect or do not connect gray matter areas in 
the dataset, respectively. We also use the Average Bundle Coverage (ABC), 
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Figure 8: Quantitative evaluation of the effect of enhancements. Evo¬ 
lution of the local error and global metrics for three different choices of £>44 and 
four different SNRs as we increase the diffusion time t. The top row shows the 
average angular error of the FOD peaks, the rows below show the average bun¬ 
dle coverage (ABC), connection to seed ratio (CSR) and the valid connection 
to connection ratio (VCCR), computed from the tractography results. 


the percentage of voxels in a bundle that is crossed by a valid streamline, av¬ 
eraged over all bundles. We combine the (VC), (1C) and (NC) in two metrics 
introduced in 


73 


• Connection to Seed Ratio (CSR), which represents the probability that 
a generated fiber actually connects two gray matter areas, computed as 
100%-NC. 


• Valid Connection to Connection Ratio (VCCR), the probability that a 
connecting fiber is correct, computed as VC/(VC-|-1C). 
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The results for the ABC, CSR and VCCR with the same enhancement pa¬ 
rameters and SNRs as for the local metric are given in Fig. Similar remarks 
hold for the global metrics as for the angular error. For all three metrics and 
all SNRs the enhancements lead to an improvement compared to CSD, the only 
exception being the ABC for SNR = 10 and D 44 = 0.02. Furthermore, as 
the SNR decreases, the larger diffusion times are beneficial and the more sig- 
nihcant the improvement is. The best results are obtained for D 44 = 0.005. 
We expect that truncation of the spherical harmonics already introduces some 
angular smoothing of the FODs on this artificial dataset, explaining the small 
effect of D 44 in the experiments. Furthermore, we see that the diffusion time t 
truly acts as a regularization parameter, resulting in a robustness for the met¬ 
rics with respect to the SNRs: the higher the diffusion time, the smaller the 
differences in the metrics between the different SNRs. 

Seeding from the white matter voxels can lead to an over-representation of 
the number of fibers in longer fiber bundles with respect to the shorter bundles 


74 . The longer bundles thereby have a larger contribution to the global metrics 


than the shorter bundles, which could lead to an overestimation of the fiber 


bundles. As proposed in 75 , we compared the global metrics when seeding from 


the gray/white matter interface for CSD and one specific set of enhancement 
parameters. The global metrics for that seeding strategy were slightly lower for 
CSD and comparable when including enhancements. For the sake of comparing 
our enhancement method with CSD, we therefore believe it is fair to use seeding 
from the white matter mask. 

The convincing improvement in the global metrics is supported by Fig. 
that shows a selection of the fiber bundles in the dataset. It can be seen that 
after enhancements, there are more valid connections in the green bundle and 
less wrong exits in the red bundle, leading to a higher (VCCR) and a better 
bundle coverage. The glyphs in the top row show that the enhancements improve 
alignment of glyphs, especially at the boundary of the fiber bundles, where the 
original CSD result tends to suffer from partial volume effects. 


3.2 Evaluation and Comparison on DW-MRI Data 

In this experiment we consider a DW-MRI dataset of a part of a human brain. 


previously used in 27 . The study was approved by the local ethical commitee 


of Maastricht University, and informed written consent was obtained from the 
subject. Although the dataset consists of only 10 axial slices, the corpus callo¬ 
sum, corona radiata and superior longitudinal fasciculus are (partly) present in 
the data. We show that the combination of CSD and enhancement is well-suited 
for different combinations of the 5-value and the number of gradient directions 
used in the acquisition. Furthermore, we make a qualitative comparison with the 


DTI-based method of 27 on this dataset and conclude with a brief quantitative 


comparison with this method on the dataset of 3.1 
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Figure 9: Comparison of CSD and enhanced CSD on the ISBI dataset. 

The top row shows that glyphs visualizing the FOD are better aligned, especially 
at the boundary of the bundle. In this case, color corresponds to the radius of 
the glyph. The bottom row depicts the tractography results, showing only the 
streamlines that pass through the indicated spheres. Here the SNR = 4 and the 
parameters used for the enhancements are £>33 = 1.0, D 44 = 0.02, t = 4.0. The 
ground truth image with the same viewpoint as the bottom figures is depicted 
on the left. 


3.2.1 Robustness with Respect to the Acquisition Parameters 

The acquisition was performed on a 3T Siemens Allegra scanner, with FOV 
208x208mm and voxel size 2x2x2mm. During the data acquisition, a brain 
region consisting of 10 axial slices was scanned with the following combinations 
of 6 -values and No, the number of orientations: b = 1000 s/mm^ with No = 
49, b = 1000 s/mm^ with No = 121 and b = 4000 s/mm^ with No = 49. 
We use again CSD with spherical harmonics up to order 8 . The higher b- 
value is obtained by using a stronger gradient pulse, making the acquisition 
more sensitive to detail in the tissue structure, but also inducing a lower SNR. 
Increasing the number of gradient directions gives a better angular resolution. 
We use deterministic tractography, with three seed regions manually selected 
in the middle of the corpus callosum, corona radiata and superior longitudinal 
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fasciculus. 

In the right column of Fig. 10 we show that after enhancements, the FOD 
allows for a more coherent reconstruction of the three bundles. Especially in 
the region where the three bundles come together, it can be seen that the fibers 
have a better propagation through the crossings. Moreover, the FODs after 
enhancements are very similar to each other, visible in the glyph visualization, 
leading to three tractography results supporting similar fiber bundles. This is 
an improvement with respect to CSD without enhancement, shown in the left 
column of Fig. |10| There we find more noisy FODs with more variation between 
the different protocols. This is also reflected in the tractography results, that 
contain more spurious fibers than after the enhancements. 

We conclude, just like in the first experiment on the phantom data, that 
applying enhancements induces more robust tractography also on real DW- 
MRI data, in this case in the sense that it is less sensitive to the acquisition 
parameters b and Ng. 


3.2.2 Comparison with a DTI-based FOD 

In the next experiment we compare the performance of our combination of CSD 


with enhancements with the method in 27 which proposed to combine DTI 


with non-linear PDE-based enhancement obtained from successively applying 
erosions and diffusions. Let us briefly describe this method, for details we 
refer to 27 , and an implementation of the PDE enhancements can be found 


in the HARDI package for Mathematica available at (http://bmia.bmt.tue. 
nl/people/RDuits/HARDIAlgorithms.zip). First an FOD on positions and 
orientations that we call Udti was constructed via a transformation of the 
tensor field D fitted to the data [^, according to the following definition 27 76 


C^DT/(y,n) = 


.(n^D-'(y)n)-T 


(23) 


dTT ^det(D(x))(is 

This FOD is then sharpened with PDE erosions, a type of morphological en¬ 
hancement adapted from 15 , on x and regularized with nonlinear diffu¬ 


sions to find crossing structures from DTI. 


Previously in 27 , the same dataset as in Fig. 10 for acquisition parameters 
b = 1000 s/mm^ and No = 49 was processed. Here we compare the FOD 
obtained with CSD, that we call Ucsd here, with Uoti in the top and bottom 
figures, respectively, of Fig. El Unlike DTI, which is limited by the Gaussian 
assumption of the diffusion profile, CSD can estimate multiple fiber orientations 
within a voxel. Furthermore, we see that the large glyphs in the Centrum 
Semiovale in the bottom figure are not apparent in Ucsd- Applying (linear) 
enhancements, as explained in Section |2(^ to Ucsd gives the second figure, and 
the approach in 27 using erosions/(nonlinear) enhancements applied to Ucti 
gives the second figure from below. It can be seen that also the enhanced DTI 
glyphs supports multiple fiber directions within voxels via extrapolation 27p8 


but at the cost of high regularization. Another noticeable difference is the fact 
that the glyphs in the CSD case are slimmer and crossings are more clearly 
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Enhanced FOD 



Figure 10: Comparison between CSD and enhanced CSD of tract- 
ography on human data. The tractography results on CSD and enhanced 
CSD data of the corpus callosum (mostly red), coronia radiata (mostly blue) 
and superior longitudinal fasciculus (mostly green), with the color related to 
the fiber direction. Enhancements are performed with D 33 = 1.0, D 44 = 0.02, 
t = 4.0. All three bundles are more apparent after enhancements and more 
fibers pass the crossings. 


defined. Whether two separate maxima are visible at a crossing is less dependent 
on the diffusion parameters in the PDE diffusion. 

Besides the visual comparison of the FOD glyphs, we provide determinis¬ 
tic tractography results for both procedures in the middle of Fig. It can 
be observed that both methods produce reasonable results, although the one 
obtained from the enhanced DTI dataset seems oversmoothed and outliers (in¬ 
dicated with the yellow arrow) can occur. This is due to the extreme diffusion 
parameters needed to perform the FOD extrapolation. We find that visually the 
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Figure 11: Comparison of glyph fields and tractography results be¬ 
tween enhanced CSD and a DTI-based FOD. Glyph visualization of an ax¬ 
ial slice of a dataset supporting the presence of the corpus callosum (mostly red), 
corona radiata (mostly blue) and superior longitudinal fasciculus (mostly green). 
Contour enhancement for CSD is performed with D 33 = 1.0, D 44 = 0.02, t = 4. 
Erosions and nonlinear diffusions for the DTI-based method are done with pa¬ 
rameters as in 27 . The tractographies corresponding to the two methods are 
shown in the middle. Outliers such as the red fiber, indicated by the arrow, 
occur due to the use of high regularization coefficients. 


combination of CSD and linear enhancements yields better tractography than 
DTI combined with erosions and nonlinear enhancements. 
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To provide a more quantitative and complete comparison of DTI, DTI and 
nonlinear enhancements, CSD and CSD with linear enhancements, we also in¬ 
clude results of the experiment in Section |3.I| for the DTI methods, see Table 
We heuristically determined good parameter settings for the nonlinear en¬ 
hancement of DTI: erosions Eq. (59)] with Du = 0.5, II 44 = 0.2, t = 2.0 
and diffusion [2^ Eq. (55)] with Du = 0.2, D 33 = 1.0, D 44 = 0.02 and t = 3. 


In Table is shown that applying enhancements for contextual regularization 
of the FOD is beneficial for both DTI and CSD. The lower the SNR, the more 
evident the improvements become. Furthermore, we see that in terms of the 
local metric, the angular error 9 of the peak orientations, the DTI methods can 
compete with the CSD based methods. However, the global metrics are signih- 
cantly higher for CSD based methods. The quantitative results on the phantom 
data in Table are in line with the qualitative comparison on real data in Fig. 

m 


Table 1: For two SNR values, the results are shown for the DTI method 
described in Section [3.2[ with or without nonlinear enhancements. We compare 
with CSD and a specific instance of enhanced CSD with parameters D 33 = 1, 
D 44 = 0.01, t = 2. For local metric 9 lower is better, for the other metrics higher 
is better. In boldface are the best results for the DTI and CSD methods. 


SNR 4 

DTI 

DTI enh 

CSD 

CSD enh 

9 (deg.) 

33.9 

15.2 

23.4 

16.3 

ABC (%) 

14.3 

18.1 

32.9 

37.9 

CSR (%) 

50.2 

54.1 

57.6 

78.2 

VCCR (%) 

17.5 

20.0 

32.9 

43.5 

SNR 10 

DTI 

DTI enh 

CSD 

CSD enh 

9 (deg.) 

23.8 

13.0 

14.9 

11.1 

ABC (%) 

15.5 

19.9 

51.6 

51.5 

CSR (%) 

69.1 

64.6 

82.8 

85.5 

VCCR (%) 

17.1 

24.3 

56.4 

57.2 


3.3 Improved Reconstruction of the Optic Radiation 


The optic radiation (OR) is a white matter fiber bundle connecting the primary 
visual cortex and the lateral geniculate nucleus (LGN), see Fig. 12 The most 


anterior part of the OR is called the Meyer’s loop (ML), of which the exact 


location is of interest for treatment of temporal lobe epilepsy 23 54 56,57 


During neurosurgery, a part of the temporal lobe is resected. To ensure that 
the OR remains intact to prevent visual field defect, it is crucial to know the 
distance from the tip of the Meyer’s loop to the Temporal Pole (ML-TP) [^, 
which shows large interpatient variability 


77 


We use DW-MRI scans of four subjects, performed on a 3.0T Philips Achieva 
MR scanner, with b = 1000 s/mm^, No = 32 and a spatial resolution of 2x2x2 
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mm. All subjects gave written informed consent; the study was approved by 
the Medical Ethics Committee of Maastricht University Medical Center (N 
43386.068). The data is acquired from healthy volunteers, and ground-truth 
ML-TP distance is not known. Therefore accuracy of this measure of our meth¬ 
ods cannot be checked, instead we focus on consistency and reproducibility. We 
apply CSD to the data to construct the FOD, with spherical harmonics up to 
order 6 requiring the estimation of 28 coefficients (as 32 directions are insuffi¬ 
cient to estimate the 45 coefficients when a spherical harmonic order 8 is used, 
when not using super-resolution as in [^). We seed from the LGN and include 
all fibers that reach the primary visual cortex. Both regions of interest are se¬ 
lected manually on a Tl-weighted image. We use probabilistic fiber tracking as 
described in Section [T3l 

We demonstrate the effect of the enhancement of CSD and the use of the FBC 
measure in Sections |3. 3. l| and |3. 3. 2[ respectively, in this relevant clinical setting. 
A quantitative comparison of the four methods CSD (O), CSD -|- enhancement 
(A), CSD -I- FBC (B) and CSD -|- enhancement -|- FBC (A-l-B) is provided in 
We show that the enhancement and/or the removal of spurious 
but in particular the combination of both methods, allows for a more 


Section |3.3.3 
fibers 


stable computation of the ML-TP distance than the original tractography result. 


3.3.1 Effect of the Enhancement of CSD on Tractography of the OR 

In this section, we apply the PDF enhancement (step A) to the CSD FOD 
as before, with parameter settings D 33 = 1, 1)44 = 0.01 and t = 2. After 
the enhancement we apply the sharpening deconvolution transform 18 and 


probabilistic tractography with 10000 streamlines. We compare the results of 
the tractography on the subjects both before and after the enhancement in Fig. 
[T^ We see that the tracking on enhanced data generally shows less spurious 
fibers, and has a better pronounced tip of the Meyer’s loop. However, the optic 
radiation is a highly curved structure, where the advantage of the enhancement 
of elongated structures cannot be fully exploited. To further reduce the spurious 
fibers, we explore our other approach in the next section. 


3.3.2 Effect of the FBC measure on Tractography of the OR 

In this section, we apply probabilistic tractography on subject 1, with 20000 
streamlines and including state of the art data scoring as in [23] (only relying 
on the data term, i.e. A = 0 in 23 Eq.(ll)]), see Fig. 14 


The kernel parameters for the coherence quantification (step B) are set to 
£>33 = 1, £>44 = 0.04 and t = 1.4 for the convolution [^. Let T be the set of 
the 1000 most anterior fibers in a tractography of the OR, that roughly form 
the Meyer’s loop. We compute the LFBC and subsequently the RFBC for all 
the fibers in T. 

Then we take emax{r) ■= maxRFBC( 7 , T), the RFBC corresponding to the 

7Gr 

“central” fiber, in the sense that it is most coherent with the fiber bundle. We 
define the filtered set T, as 
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r, :={ 7 Gr|RFBC( 7 ,r)>e}, 0 < e < ema.- (24) 

This means the parameter e acts as a threshold parameter and can be set such 
that hbers with a high spuriousness are removed. The fiber point in T^ that 
is closest to the temporal pole defines the ML-TP distance. We repeat the 
probabilistic tractography five times with the same settings on the same data, 
to qualitatively compare different stochastic realizations of the tractography 
method. The original OR reconstructions are shown in the top row of Fig. 

We observe that due to the presence of spurious hbers, the tip of the Meyer’s 
loop (indicated by the orange spheres) is estimated at different locations. When 
we set the threshold e = O.lcmaa:, removing in these cases between 6% and 8% 
of the most spurious Hbers, we obtain the results as shown in the bottom row 
of Fig. It can be seen that the resulting Hber bundles are very similar to 
each other, demonstrating less variation in the localization of the tip. 

3.3.3 Quantitative Comparisons on Four Subjects 

To support our claims of the two previous sections, we test the effect of our 
methods on the stability of the ML-TP distance under different stochastic re¬ 
alizations. Here we perform probabilistic tractography with 10000 Hbers ten 
times with the same settings, for each of the four subjects and each of the 
four methods (CSD, CSD -|- enh, CSD -|- FBC and CSD -|- enh -|- FBC). The 
FBC measure is computed from the 1000 most anterior Hbers as in the previous 
experiment and the threshold is set to e = 0.05emax- We compare the mean 
ML-TP distance and sample standard deviation determined from the tracking 
results of each of the methods. The results are summarized in the boxplots 
in Fig. |15| The Hgure strongly supports the application of the enhancements 
methods. For subjects 1-3 the ML-TP distance shows much less variation when 
including the FBC. For all subjects also (CSD -|- enh) gives more stable results 
than just CSD. Moreover, in all cases the combination (CSD -|- enh -|- FBC) 
outperforms CSD and for all but subject 1 the combined method (CSD -|- enh 
-I- FBC) also gives better results than the enhancement or FBC individually. 
It should be remarked that higher up the graph indicates a larger resection if 
used for pre-surgical evaluation, which is not necessarily positive. However, we 
prefer to have a stable and reproducible method that can be used with a safety 
margin, then a method that is more conservative, but shows large variations. 


4 Conclusions and Discussion 

We have proposed two new tools to improve alignment of Hbers in tractography 
results: (A) the combination of CSD with contextual PDF enhancements and 
(B) a Hber to bundle coherence measure to classify spurious Hbers. Both ap¬ 
proaches rely on the same contextual processing via PDFs on the space of cou¬ 
pled positions and orientations. We validate our methodology with a variety of 
experiments on synthetic and human data. 
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Figure 12: A reconstruction of the optic radiation and its positioning 
in the brain. The left figure shows how the OR is positioned in the brain, the 
close-up on the right shows how the OR wraps around the ventricular system. 
The probabilistic tractography outputs many spurious fibers. The tip of the 
Meyer’s loop, indicated by the orange sphere, is localized on a spurious hber 
and is therefore very dependent on the realization of the tractography. As a 
result, the distance from the Meyer’s loop to the Temporal pole (ML-TP) that 
is used in temporal lobe resection surgery, shows a high variation among different 
tractography outcomes. 


In the first experiment we consider a digital phantom 66 that simulates DW- 


MRI data of a challenging configuration of multiple neural-like fiber bundles for 
different noise levels, see Fig. The combination of CSD with enhancements 
and subsequent deterministic tracking was extensively tested for varying en¬ 
hancement parameters, see Fig. The enhanced FOD peaks were compared 
with the ground truth hber orientations, showing for all SNRs that the maxima 
of the enhanced FOD coincide better with the ground truth peaks than without 
application of enhancement. Also, this improvement is particularly high for very 
low SNR values. To quantitatively evaluate the impact of the enhancement on 
the tractographies we used the Tractometer evaluation system 52 . The results, 


shown in Fig. [^conhrm the beneht, for all the metrics considered, of including 
the enhancement. Also an improved stability of the metrics with respect to dif¬ 
ferent enhancement parameters is observed. Furthermore, we found that data 
with a lower SNR requires more regularization, obtained by choosing a higher 
diffusion time t in the enhancement. These quantitative evaluations of local and 
global metrics are supported by the qualitative results in Figs. [^and[^ where 
we saw that after enhancement fibers are better aligned and propagate better 
through crossings. 

The second experiment is performed on human data of a representative area 
of the brain with crossing fiber bundles. We evaluate our combination of CSD 
and enhancement for three different (single-shell) acquisition protocols, corre- 
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Figure 13: Reconstructions of the optic radiation of four subjects 
with and without use of enhancements. For all subjects, the left image 
shows the result on the original data, the right image shows the result on the 
enhanced FOD. The enhanced version generally gives less spurious fibers and 
has a more pronounced tip of the Meyer’s loop. 


spending to different ^-values and number of gradient directions. We observed, 
see Fig. that whereas tractography on CSD without enhancement showed 
notable differences between the three acquisition protocols, tractography after 
our enhancement lead to a qualitatively similar reconstruction in all cases. This 
implies that the application of enhancement in the processing pipeline makes 
the tractography results less dependent on the scanning protocol used. 

We use the same dataset and the phantom dataset to compare our method 
qualitatively and quantitatively with previous work [23[|27| in which sharpen¬ 
ing methods and nonlinear enhancement PDFs are applied to DTI. We observed 
qualitatively on real data in Fig. 11 and quantitatively in Table[^the advantage 
of CSD, that allows to use linear enhancements with less extreme regulariza¬ 
tion parameters than with the DTI based method, resulting in a more reliable 
tractography. 


30 







Figure 14: The effect of filtering spurious fibers from a probabilistic 
tractography on subject 1 in five different instances. Top row: five 
different instances of the probabilistic tractography of the OR, viewed from 
the top, selecting only the 1000 most anterior fibers. Bottom row: the result 
after filtering the most spurious fibers for each of the instances. The red sphere 
indicates the temporal pole, the white volumes represent the LGN and the 
primary visual cortex. The orange spheres are the positions with minimal ML- 
TP distance. The green sphere indicates the position of the tip averaged over 
the five tractography results, before (top) or after filtering (bottom). There is 
less variation in the position of the tip of the Meyer’s loop in the bottom row, 
i.e. after filtering, than in the top row. The fiber bundle in the left upper corner 
is the same as the one in Fig. 


12 


For our second approach to improve fiber alignment, we introduced a fiber to 
bundle coherence measure that can be used for detecting and filtering spurious 
fibers. The fiber to bundle coherence (FBC) is computed from a tractography 
based density that we constructed using the same PDE foundation as in the 
first method. As an application we considered the reconstruction of the optic 
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Figure 15: Boxplots of the ML-TP distances. For the four subjects, we 
show the mean of the ML-TP distance over ten tractography results, plus two 
standard deviations. The four different methods are indicated with different 
colors. The combination CSD -|- enh -|- FBC is the most robust in producing 
stable results. 


radiation, a fiber bundle of which the position of the anterior extent (the Meyer’s 
loop) is of interest for temporal lobe resection surgery. Accurate and stable 
localization of the tip of the Meyer’s loop is difficult due to the presence of 
spurious fibers, as shown in Fig. We demonstrated in Figs. [T^ [T4| and [TS] 
that either by enhancement of the CSD FOD, or by removing the most spurious 
fibers using the FBC measure leads to a robust probabilistic tractography. In 
particular, the combination of both methods in one pipeline allows for a more 
stable localization of the tip of the Meyer’s loop and a more stable determination 
of the Meyer’s loop to Temporal Pole distance. 

Our experiments show that our PDF enhancement methods for contextual 
processing are an effective and widely applicable tool to both enhance CSD 
data and to remove spurious fibers from tractographies. While we used CSD to 
construct an FOD, the PDF enhancement can be applied to an FOD obtained 
with any other method. We have seen that both our methods improve fiber 
alignment in tractography results and hence provide information on structural 
connectivity of the brain white matter more robustly. In the future, we aim to 
improve this framework by using data-adaptive smoothing, for example using 
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local gauge frames 78 
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